library(survey)
library(ggplot2)

## run "PCM_cleaning.R"

## Create weighted survey design object
milconf_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

#### SOCIAL DESIRABILITY: Q41 BY QUESTION ####

Q41_means <- matrix(NA, nrow = 3, ncol = 3)

Q41_means[1,1] <- svyby(~Q41A_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means[1,2] <- confint(svyby(~Q41A_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,1] * 100
Q41_means[1,3] <- confint(svyby(~Q41A_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,2] * 100

Q41_means[2,1] <- svyby(~Q41B_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means[2,2] <- confint(svyby(~Q41B_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,1] * 100
Q41_means[2,3] <- confint(svyby(~Q41B_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,2] * 100

Q41_means[3,1] <- svyby(~Q41C_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means[3,2] <- confint(svyby(~Q41C_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,1] * 100
Q41_means[3,3] <- confint(svyby(~Q41C_a, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                na.rm = TRUE))[1,2] * 100

## Plot
Q41_plot.df <- data.frame(mean = Q41_means[,1],
                          ci_lo = Q41_means[,2],
                          ci_hi = Q41_means[,3],
                          cond = c(1:3))

Q41_plot <- ggplot(Q41_plot.df) +
  geom_pointrange(aes(x = factor(cond), y = mean, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = mean, label = c("86", "94", "81")), size = 3.2,
            hjust = 1.5)
Q41_plot <- Q41_plot + xlab("") +
  scale_x_discrete(breaks = c(1:3),
                   labels = c("Q41A: \nControl Group", "Q41B: \nEveryone Supports \nthe Troops", 
                              "Q41C: \nFewer Support \nthe Troops")) +
  ylab("Confidence (% Agreement)") + ylim(50,100)
Q41_plot + theme_bw()

#### SOCIAL DESIRABILITY: Q41 BY QUESTION AND PARTY ####

Q41_means_p <- matrix(NA, nrow = 3, ncol = 9)

## Calculate means and CIs for republicans
Q41_means_p[1,1] <- svyby(~Q41A_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[1,2] <- confint(svyby(~Q41A_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[1,3] <- confint(svyby(~Q41A_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[2,1] <- svyby(~Q41B_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[2,2] <- confint(svyby(~Q41B_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[2,3] <- confint(svyby(~Q41B_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[3,1] <- svyby(~Q41C_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[3,2] <- confint(svyby(~Q41C_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[3,3] <- confint(svyby(~Q41C_a_rep, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

## Calculate means and CIs for independents
Q41_means_p[1,4] <- svyby(~Q41A_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[1,5] <- confint(svyby(~Q41A_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[1,6] <- confint(svyby(~Q41A_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[2,4] <- svyby(~Q41B_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[2,5] <- confint(svyby(~Q41B_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[2,6] <- confint(svyby(~Q41B_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[3,4] <- svyby(~Q41C_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[3,5] <- confint(svyby(~Q41C_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[3,6] <- confint(svyby(~Q41C_a_ind, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

## Calculate means and CIs for democrats
Q41_means_p[1,7] <- svyby(~Q41A_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[1,8] <- confint(svyby(~Q41A_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[1,9] <- confint(svyby(~Q41A_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[2,7] <- svyby(~Q41B_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[2,8] <- confint(svyby(~Q41B_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[2,9] <- confint(svyby(~Q41B_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

Q41_means_p[3,7] <- svyby(~Q41C_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)[1,2] * 100
Q41_means_p[3,8] <- confint(svyby(~Q41C_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,1] * 100
Q41_means_p[3,9] <- confint(svyby(~Q41C_a_dem, ~DOV_ASSIGNMENT_A, milconf_design, svymean, 
                                  na.rm = TRUE))[1,2] * 100

## Plot
Q41_means_p_plot.df <- data.frame(rep = Q41_means_p[,1],
                                  rep_lo = Q41_means_p[,2],
                                  rep_hi = Q41_means_p[,3],
                                  ind = Q41_means_p[,4],
                                  ind_lo = Q41_means_p[,5],
                                  ind_hi = Q41_means_p[,6],
                                  dem = Q41_means_p[,7],
                                  dem_lo = Q41_means_p[,8],
                                  dem_hi = Q41_means_p[,9],
                                  dcond = c(1,5,9),
                                  icond = c(2,6,10),
                                  rcond = c(3,7,11))

Q41_means_p_plot <- ggplot(Q41_means_p_plot.df) +
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2)
Q41_means_p_plot <- Q41_means_p_plot + xlab("Treatment Condition") +
  scale_x_continuous(breaks = c(2,6,10),
                     labels = c("Q41A", "Q41B", "Q41C")) +
  ylab("Confidence (% Agreement)") + ggtitle("Q41 by Question and Party ID")
Q41_means_p_plot + theme_bw()
